R Notebook

#STEP 1: Choose the variables

We aim to investigate whether COVID-19 severity influences COL11A2 expression, whether patient sex affects COL11A2 expression and whether there is an interaction between sex and disease severity in determining COL11A2 levels. We investigate three questions (H0 = null hypothesis and H1= alternative hypothesis):

H0a: The mean standardized COL11A2 expression is equal across all COVID-19 severity groups. H1a: At least one COVID-19 severity group has a different mean standardized COL11A2 expression.

H0b: The mean standardized COL11A2 expression is equal between male and female patients. H1b: The mean standardized COL11A2 expression differs between male and female patients.

H0c: There is no interaction between sex and COVID-19 severity on mean standardized COL11A2 expression. H1c: The effect of COVID-19 severity on mean standardized COL11A2 expression differs between males and females.

Why COL11A2? This protein was selected because it was mentioned in the reference article. COL11A2 is a collagen component involved in two pathways that were found to be downregulated in severe COVID-19: Integrin cell surface interactions and Collagen biosynthesis and modifying enzymes. Since these pathways include multiple collagen proteins, COL11A2 is a biologically meaningful candidate for investigating severity-associated differences.

#install.packages("ggpubr")

library(ggplot2)
library(dplyr)
library(ggpubr)
library(tidyverse)

#STEP 2: Make the dataset To test our hypotheses, we need a dataset that contains both the clinical variables (sex and WHO_temp_severity) and the standardized, log-transformed expression values of the protein of interest (COL11A2). These variables originate from different sections of the full dataset, so a new dataframe must be constructed.

During the earlier data-exploration phase, we verified that there were no missing values in sex, WHO_temp_severity or the COL11A2 protein measurements. Therefore, no additional imputation or row removal is required in this preprocessing step.

We proceed by importing the full dataset, extracting the protein expression matrix, applying a log-transformation followed by standardization (z-scaling), and combining it with the relevant clinical features into a single analysis-ready dataframe.

# Import whole dataset
combined_dataset <- read.csv("combined_dataset.csv")

# Split dataset in subsets
combined_dataset_protein <-combined_dataset[,1:7241]
proteins_log <- log(combined_dataset_protein[,2:7241])
proteins_scaled <- as.data.frame(scale(proteins_log))

feature_df <-combined_dataset[,7242:7258]

data_multivariate <- data.frame( COL11A2 = proteins_scaled$COL11A2, severity = feature_df$WHO_temp_severity, sex = feature_df$sex)
data_multivariate$severity <- factor(data_multivariate$severity, levels = c("NEGATIVE", "mild", "moderate", "severe", "critical"))

#STEP 3: Check for outliers To identify potential extreme values, we will calculate the interquartile range (IQR) of COL11A2 and determine which observations fall outside the 1.5×IQR boundaries, as these will be considered statistical outliers.

protein <- combined_dataset$COL11A2

Q1 <- quantile(protein, 0.25, na.rm = TRUE)
Q3 <- quantile(protein, 0.75, na.rm = TRUE)
IQR_value <- IQR(protein, na.rm = TRUE)

lower_bound <- Q1 - 1.5 * IQR_value
upper_bound <- Q3 + 1.5 * IQR_value

outliers <- protein[protein < lower_bound | protein > upper_bound]

outliers

 [1] 14059.1 11043.5  8828.0  8977.9 10079.4  9938.0 10522.9  8040.3  9685.6 11337.0  8650.3 11166.6 11076.8 12590.7  9094.0 11050.3  8299.4
[18] 10443.6 14506.3 13175.1  9877.9 13243.2 13882.3

ggplot(data_multivariate, aes(x = severity, y = COL11A2, fill= severity)) +
  geom_boxplot(outliers = FALSE) +
  geom_jitter (width =0.15, alpha =0.4, color = "black") +
  scale_fill_brewer(palette = "Set2") +
  theme_minimal()+
  labs(title = "COL11A2 levels across the different severity groups")

INTERPRETATION: Although several outliers are present (as seen in boxplot), there is no evidence of measurement or data-entry error and the values appear biologically plausible. Since we have no objective reason to exclude them, the outliers were retained for all subsequent analyses.

#STEP 4: Data exploration

We first visualized the distribution of log-transformed and standardized COL11A2 across the different Covid-19 severity levels. This showed that COL11A2 will be downregulated as Covid-19 becomes more severe

ggplot(data_multivariate, aes(x = severity, y = COL11A2, fill= severity)) +
  geom_boxplot(outliers = FALSE) +
  geom_jitter (width =0.15, alpha =0.4, color = "black") +
  scale_fill_brewer(palette = "Set2") +
  theme_minimal()+
  labs(title = "COL11A2 levels across the different severity groups", x = "severity groups", y = "Standardized log(COL11A2) expression")

We also looked at the distribution of COL11A2 in male and female patients. We noticed that there was no big difference between the two groups.

ggplot(data_multivariate, aes(x = sex, y = COL11A2, fill= sex)) +
  geom_boxplot(outliers = FALSE) +
  geom_jitter (width =0.15, alpha =0.4, color = "black") +
  scale_fill_brewer(palette = "Set2") +
  theme_minimal()+
  labs (title = "Sex-Based Differences in COL11A2 Expression ")

Interpretation: Only within the severe and critical groups do we observe a more pronounced difference between male and female patients. This may be explained by the limited sample size in these categories, as well as potential imbalance in sex distribution. Based on this boxplot, we could hypothesize that COL11A2 protein expression is progressively downregulated with increasing COVID-19 severity.

#STEP 5: Check assumptions

Normalization of protein COL11A2

ggqqplot(proteins_scaled$COL11A2) +
  ylim(0,5)+
  xlim(0,3)+ 
  labs(title = "Normalization of COL11A2", x= NULL, y=NULL)

The QQ-plot shows that COL11A2 values follow the theoretical normal distribution reasonably well, with only mild deviations in the upper tail.

Normalization of COL11A2 within severity groups

data_multivariate %>%
  ggplot(aes(sample = COL11A2)) +
  geom_qq() + 
  geom_qq_line() +
  facet_wrap(~severity) +
  labs( title= "Normalization of COL11A2 within severity groups" )

We assessed the normality of COL11A2 within each severity group using QQ-plots. Across groups, the distributions follow the theoretical normal line reasonably well, with mild deviations, especially in the severe and critical categories, largely attributable to small sample sizes and a few extreme observations. These deviations are not substantial enough to violate the normality assumption for the subsequent ANOVA.

Normalization of COL11A2 within sex groups

data_multivariate %>%
  ggplot(aes(sample = COL11A2)) +
  geom_qq() + 
  geom_qq_line() +
  facet_wrap(~sex)+
  labs( title= "Normalization of COL11A2 within sex groups" )

INTERPRETATION: The distribution deviates slightly from perfect normality due to the presence of a few outliers, but overall the residual pattern remains sufficiently close to normal for the assumptions of the model to hold.

Check assumptions based on a simple linear regression model

lm_COL11A2 <- lm(COL11A2 ~ sex + severity + sex*severity, data= data_multivariate)
plot(lm_COL11A2)

INTERPRETATION:

• Linearity The residuals-versus-fitted plot indicates that the linearity assumption is met. The smoothing line shows almost no curvature, and the residuals are randomly scattered around zero, suggesting no systematic under- or overprediction. The spread of residuals does not show a funnel shape, and although a few outliers are present, they do not distort the overall pattern.

• Normality of residuals The QQ-plot shows that most residuals lie close to the theoretical normal line, with only minor tail deviations. This indicates that the normality assumption is reasonably satisfied.

• Homoscedasticity (equal variance) The residual spread appears roughly constant across fitted values. The red smoothing line remains nearly flat and no funneling or curvature is observed, suggesting that the equal-variance assumption is met.

• Outliers A small number of outliers are visible, but they appear to be isolated observations rather than signs of structural problems in the data. Because no evidence of measurement error was found, they are retained in the analysis.

STEP 6: The two-sample Anova test

Because our data is approximately normally distributed, we selected a two-way ANOVA to test our multivariate hypothesis.

anova_model <- aov( COL11A2 ~ severity + sex + severity * sex , data = data_multivariate)
summary(anova_model)

              Df Sum Sq Mean Sq F value Pr(>F)    
severity       4  89.52  22.380  34.804 <2e-16 ***
sex            1   0.70   0.705   1.096  0.296    
severity:sex   4   1.23   0.309   0.480  0.751    
Residuals    231 148.54   0.643                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

library(effectsize)

eta_squared(anova_model)

# Effect Size for ANOVA (Type I)

Parameter    | Eta2 (partial) |       95% CI
--------------------------------------------
severity     |           0.38 | [0.29, 1.00]
sex          |       4.72e-03 | [0.00, 1.00]
severity:sex |       8.24e-03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

INTERPRETATION: At the 5% significance level, the two-way ANOVA shows a statistically significant effect of COVID-19 severity on standardized log-transformed COL11A2 expression (F(4, 231) = 34.8, p < 0.001). In contrast, neither sex (F(1, 231) = 1.10, p = 0.296) nor the interaction between sex and severity (F(4, 231) = 0.48, p = 0.751) is significant. We therefore reject the null hypothesis for severity, but do not reject the null hypotheses for sex and for the sex-severity interaction.

effect size: COVID-19 severity shows a large effect on COL11A2 expression (partial η² = 0.38), indicating that approximately 38% of the variance is explained by severity. In contrast, sex and the sex×severity interaction have negligible effect sizes (partial η² < 0.01).

The anova model was repeated for the interaction between COL11A2 and severity to use for the post hoc analysis

sign_anova <- aov( COL11A2 ~ severity , data = data_multivariate)
summary(sign_anova)

             Df Sum Sq Mean Sq F value Pr(>F)    
severity      4  89.52  22.380    35.1 <2e-16 ***
Residuals   236 150.48   0.638                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#STEP 7: Post hoc analyse

We know that Covid-19 severity has a significant effect on the expression of COL11A2. It is important that we do a Post-hoc test to find out for which groups the expression is significantily different.

#install.packages("multcomp")
library(multcomp)

mcp <- glht(sign_anova,linfct = mcp (severity = 'Tukey' ))
summary(mcp)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = COL11A2 ~ severity, data = data_multivariate)

Linear Hypotheses:
                         Estimate Std. Error t value Pr(>|t|)    
mild - NEGATIVE == 0      -0.5923     0.1354  -4.373 0.000168 ***
moderate - NEGATIVE == 0  -1.1851     0.1583  -7.489  < 1e-04 ***
severe - NEGATIVE == 0    -1.7688     0.1937  -9.130  < 1e-04 ***
critical - NEGATIVE == 0  -2.0581     0.2471  -8.328  < 1e-04 ***
moderate - mild == 0      -0.5928     0.1390  -4.264 0.000249 ***
severe - mild == 0        -1.1764     0.1784  -6.595  < 1e-04 ***
critical - mild == 0      -1.4658     0.2353  -6.230  < 1e-04 ***
severe - moderate == 0    -0.5836     0.1963  -2.974 0.024695 *  
critical - moderate == 0  -0.8730     0.2491  -3.504 0.004543 ** 
critical - severe == 0    -0.2893     0.2730  -1.060 0.817948    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

confint(mcp)$p

NULL

INTERPRETATION: The Tukey post-hoc test shows that mean standardized log-COL11A2 expression is significantly lower in all COVID-19 severity groups (mild, moderate, severe, critical) compared with the NEGATIVE group (all adjusted p < 0.001). For example, the critical group has on average 2.06 standard units lower COL11A2 expression than the NEGATIVE group (estimate = −2.06, 95% CI [−2.73, −1.39]). In addition, each successive increase in severity from mild to moderate, severe and critical is associated with a significant further decrease in COL11A2 expression (all adjusted p < 0.05), except for the contrast between severe and critical, where the 95% confidence interval includes zero (estimate = −0.29, 95% CI [−1.03, 0.45], p = 0.82). At a global 5% significance level we therefore conclude that COL11A2 expression decreases progressively with increasing COVID-19 severity.

Visualisation

This boxplots supports the results obtained via the two-way ANOVA test

library(ggplot2)


ggplot(data_multivariate, aes(x = severity, y = COL11A2, fill = sex)) +
  geom_boxplot(outliers = FALSE) +
  geom_jitter (width =0.15, alpha =0.4, color = "black") +
  scale_fill_brewer(palette = "Set2") +
  theme_minimal() +
  labs(title = "The effect of sex and Covid-19 severity on COL11A2 expression", x = "severity groups", y = "Standardized log(COL11A2) expression")

FINAL CONCLUSION:

H0a (severity): We reject the null hypothesis (p < 0.001). COVID-19 severity has a significant effect on standardized log-COL11A2 expression. Both the ANOVA and Tukey post-hoc tests show a progressive decrease in COL11A2 levels from NEGATIVE → mild → moderate → severe → critical, with all pairwise differences significant except severe–critical. The boxplot visually confirms this downward trend, indicating that COL11A2 is strongly downregulated as disease severity increases.

H0b (sex): We do not reject the null hypothesis (p = 0.30). There is no statistically significant difference in mean COL11A2 expression between male and female patients. The figure shows substantial overlap between sexes in each severity category, supporting this result.

H0c (interaction): We do not reject the null hypothesis (p = 0.75). The effect of COVID-19 severity on COL11A2 expression does not differ between males and females. The boxplot shows small apparent sex differences only in the severe and critical groups, but these patterns are likely due to limited and imbalanced sample sizes rather than a true interaction.

Overall conclusion: Overall, COL11A2 expression decreases consistently with increasing COVID-19 severity, indicating a strong and significant severity effect visible in both the statistical tests and the final figure. Sex and the sex×severity interaction are not significant, and the minor visual differences in the most severe groups likely reflect small sample sizes rather than true biological effects. Thus, the severity-related decline in COL11A2 appears robust and independent of patient sex.